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Abstract. In this paper, we propose a geometric integrator for nonholonomic 
mechanical systems. It can be apphed to discrete Lagrangian systems specified 
through a discrete Lagrangian : Q X Q — > M, where Q is the configuration 
manifold, and a (generally nonintegrable) distribution D C TQ. In the pro- 
posed method, a discretization of the constraints is not required. We show 
that the method preserves the discrete nonholonomic momentum map, and 
also that the nonholonomic constraints are preserved in average. We study 
in particular the case where Q has a Lie group structure and the discrete La- 
grangian and/or nonholonomic constraints have various invariance properties, 
and show that the method is also energy-preserving in some important cases. 
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1. Introduction 

During the last years, there has been an increasing interest in nonholonomic 
mechanical systems, in part motivated by some open questions in the subject, such 
as those concerning reduction, integrability, stabilization or controllability; and also 
for their applicability in engineering, specially in robotics, mainly since it describes 
the motion of wheeled devices (see [31[31[TT] and the expository paper jlj). 

When a mechanical system is subjected to some external constraints, the lat- 
ter may be expressed in terms of relations imposing restrictions on the allowable 
positions and velocities. The constraints are then called nonholonomic if the ve- 
locity dependence is essential, in the sense that the constraint relations can not be 
reduced, by integration, to relations depending on the position coordinates only. 
Geometrically, nonholonomic constraints are globally described by a submanifold 
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V of the velocity phase space TQ. In most of the known examples I? is a vector 
subbundle of TQ, i.e., the constraints have a linear dependence on the velocities. 
Lagrange-d'Alembert's principle allow us to determine the set of possible values 
of the constraint forces from the constraint manifold T). Then, to determine the 
dynamics of the nonholonomic system, it is only necessary to fix initially the pair 
(L,I?), where L : TQ ^ K is a Lagrangian function, usually of mechanical type (see 
[2l|6l|8] for an extension of the classical Lagrange-d'Alembert's principle). 

Very recently, many authors [lOl [121 1131 HSl [22] started the study of geomet- 
ric integrators adapted to nonholonomic systems, obtaining very stable numerical 
integrators with some preservation properties (such as discrete nonholonomic mo- 
mentum map preservation) and very good energy behavior. This problem is of 
considerable interest given the crucial role of nonholonomic dynamics in many ap- 
plications in engineering. From the numerical point of view, in |23j it appeared 
as an open question: "...The problem for the more general class of non-holonomic 
constraints is still open, as is the question of the correct analogue of symplectic 
integration for non-holonomically constrained Lagrangian systems...". 

The most interesting approach to nonholonomic integrators appears as an adap- 
tation of the so-called variational integrators 2\\ incorporating a discrete constraint 
submanifold, in addition to a discretization of the Lagrangian function and the 
vector subbundle T>. Then, the numerical method is obtained from the so-called 
Discrete Lagrange-d'Alembert's principle [12], recovering many of the geometric 
properties of the continuous system. 

Obviously, since nonholonomic mechanics is not symplectic-preserving, it seems 
interesting to try to preserve another geometric invariance property of the contin- 
uous nonholonomic system, as for instance, the energy function in the autonomous 
case. This is precisely the starting point of view of our paper. Moreover, a dis- 
cretization of the constraints is not required here. We show that the method pre- 
serves the discrete nonholonomic momentum map, and also that the constraints are 
preserved in average. We study in particular the case where the configuration space 
is a Lie group and the discrete Lagrangian and/or nonholonomic constraints have 
various invariance properties, and show that the method is also energy-preserving 
in many important cases. In particular, the main result of the paper. Theorem [T] 
states that if the configuration space is a Lie group and the Lagrangian is defined 
by a bi-invariant Riemannian metric, then, from a left-invariant discretization of 
the Lagrangian, we obtain a fixed time-step, energy-preserving numerical 
method for the continuous nonholonomic system, without requiring any invari- 
ance conditions on T>. See [9 for a variable time-step algorithm that preserves 
energy. 

The paper is structured as follows. In Section [2] we introduce continuous non- 
holonomic mechanical systems for the case of mechanical energy Lagrangians de- 
fined by a given Riemannian metric and a potential function. In this case, the 
equations of motion for the constrained system are geodesic equations for an afline 
connection (in the kinetic case) that is not generally Levi-Civita, obtained from 
the induced orthogonal projection onto the nonholonomic distribution (see [TlllTj'l. 
In Section [Sj we recall some definitions concerning discrete variational mechanics 
(discrete Lagrangian, discrete Euler-Lagrange equations, discrete flow, momentum 
map...). The new proposed method appears in Section[4] constructed from the dis- 
crete Lagrangian and the orthogonal projectors induced by the distribution T> and 
the Riemannian metric. Then we consider the case when the configuration space is 
a Lie group and we obtain under adequate invariance properties the preservation of 
energy. In addition, we study the momentum nonholonomic map for the proposed 
nonholonomic integrator. In Section [5| we introduce a nonholonomic version of the 
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Stormer-Veiiet method which is a natural extension of the RATTLE method for 
nonholonomic systems. In Section [6] we test our method in three examples (the 
nonholonomic particle, the snakeboard and the Chaplygin sleigh). The paper ends 
with a section of conclusions and future work. 

2. Continuous nonholonomic mechanics 

We shall start with a configuration space Q, which is an n-dimensional difFeren- 
tiable manifold with local coordinates (g'), 1 <i < n = dimQ. Constraints linear 
in the velocities are given by equations of the form 

r{q\ql^t^nqW^O, l<a<m, 

depending, in general, on configuration coordinates and their velocities. From an 
intrinsic point of view, the linear constraints are defined by a distribution V on Q 
of rank n — m such that the annihilator of V is locally given by 

V° = span {/x" = fi°;dq' ; 1 < a < m} 

where the one- forms /i° are independent. 

The various kinds of constraints we are concerned with will roughly come in 
two types: holonomic and nonholonomic, depending on whether the constraint 
is derived from a constraint in the configuration space or not. Therefore, the 
dimension of the space of configurations is reduced by holonomic constraints but 
not by nonholonomic constraints. Thus, holonomic constraints allow a reduction in 
the number of coordinates of the configuration space needed to formulate a given 
problem (see j^). 

We will restrict ourselves to the case of nonholonomic constraints. In this case, 
the constraints are given by a nonintegrable distribution T). In addition to these 
constraints, we need to specify the dynamical evolution of the system, usually by 
fixing a Lagrangian function L : TQ — > K. In mechanics, the central concepts 
permitting the extension of mechanics from the Newtonian point of view to the 
Lagrangian one are the notions of virtual displacements and virtual work; these 
concepts were formulated in the developments of mechanics, in their application 
to statics. In nonholonomic dynamics, the procedure is given by the Lagrange- 
d'Alembert principle. This principle allows us to determine the set of possible 
values of the constraint forces from the set V of admissible kinematic states alone. 
The resulting equations of motion are 

-d (dL\ ^L^ 

where Sq^ denotes the virtual displacements verifying 

V - 

(for the sake of simplicity, we will assume that the system is not subject to non- 
conservative forces). This must be supplemented by the constraint equations. By 
using the Lagrange multiplier rule, we obtain 

dt\dq^J dt^ 

The term on the right represents the constraint force or reaction force induced 
by the constraints. The functions Ao are Lagrange multipliers which, after being 
computed using the constraint equations, allow us to obtain a set of second order 
differential equations. 
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Now we restrict ourselves to the case of nonholonomic mechanical systems where 
the Lagrangian is of mechanical type 

1 

T 

where <? is a Riemannian metric on the configuration space Q. Locally, the metric 
is determined by the matrix M = {gij)i<ij<n where gij = g{d/dq^,d/dq^). 

Using some basic tools of Riemannian geometry, we may write the equations of 
motion of the unconstrained system as 

Ve(t)c(t) = -grad V{c{t)), (1) 

where V is the Levi-Civita connection associated to g. Observe that if y = 
then the Euler-Lagrangian equations are the equations of the geodesies for the 

Levi-Civita connection. 

When the system is subjected to nonholonomic constraints, the equations become 

Vc(t)c(t) = -grad V{c{t)) + X{t), c{t) £ 

where A is a section of V-^ along c. Here V-^ stands for the orthogonal complement 
of V with respect to the metric g. 

In coordinates, by defining the functions T^j (Christoffel symbols for V) by 

we may rewrite the nonholonomic equations of motion as 

q'^{t) + r':^{c{tM{tW{t) = -g'\c{t))—^ + Xa{t)g'Mt))f,t{c{t)) 

where t ^ {q^{t),. . . , q'"{t)) is the local representative of c and {g^^) is the inverse 
matrix of M. 

Since 5 is a Riemannian metric, the mxm matrix (C"**) = {jJ-i g^-^ iJ^'j) is symmetric 
and regular. Define now the vector fields Z", 1 < a < m on Q by 

giZ", Y) = for all vector fields Y, 1 < a < m; 

that is, Z° is the gradient vector field of the 1-form /x". Thus, T)-^ is spanned by 
Z"", 1 < a < m. In local coordinates, we have 

d 

^ dq^ 

We can construct two complementary projectors 

orthogonal with respect to the metric g. The projector Q is locally described by 

Q = CabZ'^ ® = Cabg^'fiffil^ ® dg^ 

Using these projectors we may rewrite the equations of motion as follows. A curve 
c{t) is a motion for the nonholonomic system if it satisfies the constraints, i.e., 
c{t) e ■E'c(t)) and, in addition, the "projected equation of motion" 

^(Vc(t)c(i)) = -P(grad V{cm (2) 

is fulfilled. 
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Summarizing, we have obtained the dynamics of the nonholonomic system ^ 
applying the projector V to the dynamics of the free system ([T]). In Section |4] we 
will use V and Q to obtain a geometric integrator for nonholonomic systems. 

3. Variational integrators 

The equations of motion for an unconstrained Lagrangian system given by a 
Lagrangian function L : TQ M are the well-known Euler-Lagrange equations 

d fdL\ dL ^ , 

- TT-^ = 0, \ <i<n. 



dt \dq 

It is well known that the origin of these equations is variational (see Now, 
variational integrators retain this variational character and also some of the geo- 
metric properties of the continuous system, such as symplecticity and momentum 
conservation (see [TH HI] and references therein) . 

In the following we will summarize the main features of this type of numerical 
integrators. A discrete Lagrangian is a map L^: Q ^ Q — > M, which may be 
considered as an approximation of a continuous Lagrangian L : TQ — > M. Define 
the action sum Sd ■ Q^^^ K corresponding to the Lagrangian Ld by 

N 

Sd = y^-Ld(gfc-i,gfc), 
fe=i 

where Qk ^ Q for < fc < A^. The discrete variational principle states that the 
solutions of the discrete system determined by Ld must extremize the action sum 
given fixed endpoints go E^nd qn- By extremizing Sd over qk, 1 < k < N — 1, we 
obtain the system of difference equations 

DiLd{qk,qk+i) + D2Ld{qk-i,qk) = 0. (3) 

or, in coordinates, 

^(9fc,gfe+i) + |^('7fe-i,'Zfc) = 0, l<i<n, l<k<N-l. 
dql dq\ 

These equations are usually called the discrete Euler Lagrange equations. 
Under some regularity hypotheses (the matrix (£'12^^(9/01 9fe+i)) is regular), it is 
possible to define a (local) discrete flow TiQxQ^QxQ, by T{qk-i,qk) = 
{qk,qk+i) from Define the discrete Legendre transformations associated to Ld 
as 

W-Ld-. QxQ^T*Q 

(qo-.qi) I — > {qo,-DiLd{qo,qi)) 
V+Ld-. QxQ^T*Q 

{qo,qi) I — * {qi,D2Ld{qo,qi)) , 

and the discrete Poincare-Cartan 2-form ojd = {V~^Ld)*ujQ = {¥~ Ld)*u!Q, where 
ojQ is the canonical symplectic form on T*Q. The discrete algorithm determined 
by T preserves the symplectic form Ud, i.e., T*^^^; = ujd- Moreover, if the discrete 
Lagrangian is invariant under the diagonal action of a Lie group G, then the discrete 
momentum map Jd : Q x Q — > g* defined by 

{Jd{qk,qk+i),0 = {D2Ld{qk,qk+i)AQ{qk+i)) 

is preserved by the discrete flow. Therefore, these integrators are symplectic- 
momentum preserving. Here, denotes the fundamental vector field determined 
by ^ Eg, where g is the Lie algebra of G. 
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4. A GEOMETRIC NONHOLONOMIC INTEGRATOR 

This work proposes a numerical method for the integration of nonholonomic 
systems. It is not truly variational; however, it is geometric in nature and we show 
in Corollary [4] that it preserves the discrete nonholonomic momentum map in the 
presence of horizontal symmetries. Moreover, we prove in Theorem [l] that under 
certain symmetry conditions, the energy of the system is preserved. 

Consider a discrete Lagrangian Ld'- Q x Q ^ M. The proposed discrete nonholo- 
nomic equations are 

Vl{DM<lk,qk+i)) + Vl{D2Ld{qk^^,qk)) - (4a) 

Ql^{DrLd{qk,qk+l)) - Ql{DMqk^,,qk)) = 0, (4b) 

where the subscript emphasizes the fact that the projections take place in the 
fiber over qk- The first equation is the projection of the discrete Euler-Lagrange 
equations to the constraint distribution V, while the second one can be interpreted 
as an elastic impact of the system against D (see [H]). This is what will provide 
the preservation of energy. Note that we can combine both equations into 

DiLdiqk,qk+i) + {V* - Q*)i?2id((Zfc-i, %) = 0, 

from which we see that the system defines a unique discrete evolution operator if 
and only if the matrix {Di2Ld) is regular, that is, if the discrete Lagrangian is 
regular. Locally, the method can be written as 

DiLd{qk,qk+i) + D2Ld{qk-i,qk) = {'>^k)blJL^ (5a) 

9'\qk)^^nqk) (^{qk,qk+i) - ^{qk-i,qk)] - 0. (5b) 
V 9qi dq{ J 

Using the discrete Legendre transformations defined above, define the pre- and 
post-momenta, which are covectors at q^, by 

Pfe_i,fc =p+((7fe-i,(7fe) = ¥+Ld{qk-i,qk) = D2Ld{qk~i, qk) 

Pk.k+i =p^{qk,qk+i) = ^^Ld{qk,qk+i) = -DiLd{qk,qk+i)- 



In these terms, equation (5b I can be rewritten as 

g^qkjf^iiqk) I ^ I = 

which means that the average of post- and pre-momenta satisfies the constraints. 
In this sense the proposed numerical method also preserves the nonholonomic con- 
straints. 

We may rewrite the discrete nonholonomic equations as 

Plk+i^CP-QYlipt-i^k)- (6) 

We interpret this equation as a jump of momenta during the nonholonomic evo- 
lution. Compare this with the condition pj^ ^.^i — pt-i k imposed by the discrete 
Euler-Lagrange equations (that is, for unconstrained systems). In our method, 
the momenta are related by a reflection with respect to the ima ge of the projector 
V* : T*Q (X»-L)°. This is illustrated, in the context of Section [II] in figure [l] 



4.1. Left-invariant discrete Lagrangians on Lie groups. Consider a discrete 
nonholonomic Lagrangian system on a Lie group G, with a discrete Lagrangian 
Xrf : G X G ^ M that is invariant with respect to the left diagonal action of G 
on G X G (see [5J HH])- We do not impose yet any invariance conditions on the 
distribution V. If we write Wk = g^^gk+i, then we can define the reduced discrete 
Lagrangian /<j: G M as ld{Wk) = Ld{gk,9k+i)- Note that Dld{Wk) G T^fi- 
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Figure 1. Evolution of momenta, depicted here as sohd arrows. 
The right translations are a consequence of the left-invariance of 
Ld, and the reflection at gk is the proposed method. 

Computing the derivative, we obtain 

Pk,k+i ^ -DiLd{gk,9k+i) = L*^-iRlY^Dld{Wk), 

where L* and R* are the mappings on T* G induced by left and right multiplication 
on the group, respectively (this should not be confused with the Lagrangian L). 
We use this to write 

Pt.k+i = D2Ld{gk,9k+i) = L* iDldiWk) = L* iR*^-iLl^plj^^^ = R*w-'Pk,k+i- 

k k k k 

Therefore, the discrete nonholonomic equations ([6| become 

Pk.k+i -CP- QT {Rw-^Pk-i,k) ■ (7) 

The relationships between the pre- and post-momenta are depicted in figure [l] 

Note that we do not need here that the metric used to build the projectors is 
the metric giving the kinetic energy in the Lagrangian. 

4.2. Left-invariant Lagrangian and projectors. Take a left-invariant discrete 
Lagrangian L^: G x G ^ M as in the previous section, and assume that 1) and 
T)-^ are left-invariant. This is typically a consequence of V and the metric on G 
being left-invariant, although it can be assured by weaker conditions on the metric 
(preserving the orthogonality of and by left translations) . This is equivalent 
to the left-invariance of the projectors V and Q, which in turn is equivalent to the 
left-invariance of V ~ Q, as a straightforward verification shows. 

Since our goal is to rewrite equation ([7| on the dual g* of the Lie algebra, we 
define the discrete body momentum pk'. G x G q* as 

P^ = ^IkPk.k+i^ 
which agrees with the definition in [13 . Then Q reads 

Since {V — Q)* is left-invariant, we obtain 



that is, 



Pk = {V -QY {P^d*w,_,Pk-i 



4.3. Preserving energy on Lie groups. Let us now consider the case where Q 
is a Lie group G, the nonholonomic distribution T) is not necessarily G-invariant, 
and L is regular and bi-invariant. 

Since we are restricting ourselves to Lagrangians of mechanical type, the poten- 
tial energy is necessarily zero. The left-invariance of L implies that it must be of 
the form 

L{vg)^\{lg-^vg,g-\g), (8) 
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where I : g — > g* is a symmetric non-singular inertia tensoiij The bi-invariance, 
however, imposes the equivariance condition Ad*-i o I = I o Adg for all g G G, 
as is straightforward to check. We remark that in this section, the metric used 
to build the projectors will be the same that defines the Lagrangian. If we take a 
discretization : G x G M. (which needs to be fe/i-invariant only) , the equations 
of motion ([t]) hold. Then we can prove the following result. 

Theorem 1. Consider a nonholonomic system on a Lie group with a regular, hi- 
invariant Lagrangian and with an arbitrary distribution V, and take a discrete La- 
grangian that is left-invariant. Then the proposed discrete nonholonomic method Q 
is energy-preserving . 

Proof. The equivariant inertia tensor I induces an Ad-invariant scalar product on 
Q and a bi-invariant metric on G. It also defines an inner product and a 
corresponding norm || • \\i on each fiber of T*G that inherit this bi-invariance. If 
p I— > is the index-raising operation associated to the kinetic energy metric, then 

ibfflll = {Pg^PD ^ {Pg^Lgr^L*gPg) = {Pg, Rgl'^ R*gPg) . 

The dual applications of the projectors V and Q are orthogonal complementary 
projectors with respect to this inner product, and thus for p gT*G, 

\\{V-Qrp\\l^{V*p,V*p)i + {Q*p, Q*p)i^\\{V+Qrp\\l^\\p\\l 

The energy function is given in the continuous setting hy H ^ {dL/dg,g) — L 
as a function of the position g and momentum p = dL/dg. For L given by (|8| we 
have 

Hig,p)^^-{L;p,I-'L;p)^^-\\p\\l 

Proving that the energy is preserved amounts to showing that equation ([7| preserves 
II • III. Since || • ||i is in particular right-invariant, then ^J^-i '■ Tg^__^G — s- Tg^G is an 

isometry. In addition, we have shown above that (7^ — Q)* is also norm- preserving, 
so we obtain 

H{9k,Pk.k+i) = Higk-i.p'k^i^k)- n 

Remark 2. While the proof above shows that the norm of the post-momenta is 
preserved, the norm of the pre-momenta is also preserved since they are related by 
a reflection (equation ([6])). 

4.4. The average momentum. Take a discrete nonholonomic system on G as 
in the previous section, but add the condition that V is right-invariant. Since the 
metric on the group is right-invariant, so is the projector V . Take a trajectory of 
the system and define at each t/j, the average momentum 

Pfe = ^ (Pfc-i,fc+Pfc,fc+i) • (9) 
Using ([6]), Q and the fact that {V — Q)* is its own inverse, we have 

Pu = ^ ((^ - QnPu.k+i)+PKk+i) = T*iPlk+i) - T'*iR*w-2Pk-i,k) 
= R*w-iJ^*^Pk-iM) = R*w-iPk-i- 

Since the norm || • ||i on each fiber of T*G defined in the proof of Theorem [l] is 
right- invariant, we obtain |jpA;|ji — \\pk-i\\i, so 

H{gk,Pk) = H{gk-i,Pk^i)- 
In addition, by equation we have that Q*{pk) = 0, sopk satisfies the constraints. 



In the context of Lie groups, g will denote an element of G instead of the metric. 
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4.5. Preservation of the nonholonomic momentum map. Let us recall some 
concepts regarding symmetries of nonholonomic systems. Suppose that a Lie group 
G acts on the configuration manifold Q. Define, for each q G Q, the vector subspace 
q'^ consisting of those elements of g whose infinitesimal generators at q satisfy the 
nonholonomic constraints, i.e.. 

The (generalized) bundle over Q whose fiber at q is is denoted by g^. 

A horizontal symmetry is an element ^ G such that Cq('?) ^ for all q €z Q. 
Note that a horizontal symmetry is related naturally to a constant section of g^. 

Now consider a discrete Lagrangian L^,: Q x Q ^ R, and define the discrete 
nonholonomic momentum map J2^ ■ Q x Q ^ (fl^)* in [10] by 

Jf g'^ 

£, i-^ {D2La{qk-i,qk),£,Q{qk)) ■ 

For any smooth section ^ of g^ we have a function (J^^)^- Q x Q ^ R, defined as 
{Jd'^)liQk~i,qk) = Jf\qk-i,qk) (Jiqk)j ■ We can now prove the following result. 

Theorem 3. Assume that Ld is G-invariant, and let £^ be a smooth section of . 
Then, under the proposed nonholonomic integrator, (Jj'')^ evolves according to the 
equation 

{Jf)^{qk, qk+l) - (Jf gfc) - {D2Ldiqk, qk+l), " ^k)Q (qk+l)) 

where G g are the result of dropping the base points of ^(qk) and C('?fe+i) 

respectively. 

Proof. By the invariance of Ld we have 

Ld{exp{s^k)qk,exp{s^k)qk+i) = Ld{qk,qk+i), 
and differentiating at s = we get 

{DiLdiqk, qk+i) , {^k)Q{qk)) + {D2Ld{qk, qk+l) , {^k)Q{qk+i)) = 0. 
On the other hand, the proposed integrator implies 

{V ~ Qy{DiLd{qk,qk+l)) + D2Ld{qk-l,qk) ^ 0- 
From this, and using the fact that (Cfc)Q('?fe) G 2?, we have 
{Jd'^)^{(lk-i,qk) = {D2Ld{qk-i,qk),{£,k)Q{qk)) 

= - {DiLd{qk,qk+i), {V - Q) {{£,k)Q{qk))) = - {DiLd{qk,qk+i),{£.k)Q{qk)) 

= {D2Ld{qk,qk+i),{^k)Q{qk+i)) ■ 

Then 

iJf)^iqk,qk+l) - {jTy^{qk-l,qk) = 

= {D2Ld{qk,qk+i),{ik+i)Q{qk+i)) - {D2Ld{qk,qk+i)A^k)Q{qk+i)) 
= {D2Ld{qk,qk+i), {S.k+1 - S.k)Q{qk+i)) ■ □ 

Corollary 4. // Ld is G-invariant and is a horizontal symmetry, then the pro- 
posed nonholonomic integrator preserves {J^^^)^. 
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5. A THEORETICAL EXAMPLE: NONHOLONOMIC VERSION OF THE 
StORMER-VERLET METHOD 

Consider a continuous nonholonomic system determined by the mechanical La- 
grangian L : M^" ^ 

L{q,q)^]^fMq~V{q) 

(with M a constant, invertible matrix) and the constraints determined by iJ-{q)q = 
where fi{q) is a m x n matrix with rank ii — m. 
Consider now the symmetric discretization 



Ld{qk,qk+ij = ^hL Iqk, 



^ ,r f *+l - Qk 



^ {qk+i - qkf M {qk+i - gfc) - ^ {Vi^k) + V{qk+i)) 



After some straightforward computations we obtain that equations (5a I and (5b) 
for the proposed nonholonomic discrete system are 



qk+i - 2qk + qk-i = -h M [Vq{qk) + A* (gfe)Afe 

'?fc+i - Qk-i 



= Kqk) 



2h 



(lOa) 
(10b) 



where Vq{q) — {dV / dq^{q)) and the Lagrange multipliers relate to those in equation 
(5a I by = Xk/h. We recognize this set of equations as an obvious extension of 



the SHAKE method proposed by [25] to the case of nonholonomic constraints. 
The SHAKE method is a generalization of the classical Stormer-Verlet method in 



presence of holonomic constraints. Equations ( 10 1 were proposed by R. McLachlan 
and M. Peiimutter [52] (see equations (5.3) therein) as a reversible method for 
nonholonomic systems not based in the Discrete Lagrange-d'Alembert principle. 

The momentum components are approximated by the average momentum pk = 
M{qk+i - qk-i)/2h given by equation ([9|. Denoting Pk+1/2 = M{qk+i - qk)/h, 
equations (10a) and (10b I are now rewritten in the form 



Pk+1/2 =Pk-^ (Vg((7fe) + fJ^'^{qk)Xk^ , 

Qk+i = qk + 

= fi{qk)M^'^pk. 



The definition of Pk+i requires the knowledge of qk+2 and, therefore, it is is 



natural to apply another step of the algorithm ( 5a I and ( 5b I to avoid this difficulty. 
Then, we obtain the new equations: 



Pk+l = Pk+1/2 - ^ (Vq{qk+i) + li^{qk+l)\k 
= ii{qk+i)M^^pk+i. 



+1 ' 



The interesting result is that we obtain a natural extension of the RATTLE 
algorithm for holonomic systems to the case of nonholonomic systems. Unifying 
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the equations above we obtain the following numerical scheme 



2 i^qilk) + fJ- iqk)Xk 



Pk+l/2 = Pk 

Qk+i = Qk + Vfc+i/2 
h 



Pk+l — Pk+l/2 



Vq{qk+i) + ^^ ((7fc+i)Afe^ 



= /i((7fe+i)M Vfc+i- 



(11a) 
(lib) 
(He) 

(Ud) 
(lie) 



These equations allow us to take a triple {qk,Pk, Afe) satisfying the constraint equa- 
tions (11c), compute Pk+1/2 using ( |lla[ ) and then qk+i using (lib). Then, equa- 
tions (lid I and (lie) are used to compute the remaining components of the triple 
{qk+i,Pk+i,^k+i)- Of course, from Theorem[l]we obtain that, in the case V = 
the numerical method is energy preserving. 

Remark 5. From this Hamiltonian point of view, we have shown that the initial 
conditions for this numerical scheme are constrained in a natural way ((goiPo) with 
IJL{qo)M^^Po = 0), that is, the initial conditions are exactly the same as those for the 
continuous system. However, if we want to maintain the algorithm in the cartesian 
product Q X Q, then the appropriate set of initial conditions is now 

{(go, 91) e Q X Q I ¥-Ldiqo,qi) e (1?^)°} 



Mo 



(12) 



(90, 91) (^Qy-Q 



9''iqo)pnqo)^iqo,qi) 







In the particular case of the nonholonomic projection of the Stormer-Verlet method 
we have that 



A4o = |(9o,gi)GK"xIR" 
Thus, if (90,91) G Mo, we define 



li{qo)M- 



M 



91 - 90 



h 



n(9o; 







Pa 



dLd, , 91 
^—(90,91) = — 



— + ^^9(90) = Po+1/2 + ^Vq(go) 



From the expression of Mo we have that (11c I holds for fc = 0, and the definition 
of Po+1/2 yields precisely equation (lib). If we take A = then (Hal holds too. 



Thus, (9o,Po,0) can be used to initialize the algorithm (111 



Remark 6. In the particular case where the constraints are integrable, that is, the 
motion is only defined on a submanifold N of Q, then the most natural choice is 
to restrict the discrete Lagrangian to x A^: {Ld)\NxN (see [51] and references 
therein). In a local description N is determined by the vanishing of a family of 
independent functions .9^(9) — 0, 1 < a < m. Differentiating, we obtain new 
constraints 

^(9)9' = (13) 

which are satisfied by the trajectories {c{t),c(t)) in the continuous problem. 

If we directly apply our method to a holonomic system we obtain the preser- 



vation of constraints (13) but the computed numerical solution will not usually 
lie on the constraint submanifold 5^(9) = 0. For instance, it seems more natu- 
ral to change (11c I by g°'{qk+i) = 0, as appears in the classical RATTLE method. 
Nevertheless, in the case V — 0, our method has as an additional feature the preser- 
vation of energy. We could say that the proposed method is specifically designed 
for nonintegrable constraints. 
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6. Numerical examples 



Example 1. The following typical example will illustrate some of the constructions 
of previous sections. It corresponds to a discretization of the nonholonomic particle 
in described by 

L{x, y, z, X, i) = i (±2 + + i;2^ 

and the nonholonomic constraint ip ~ z — yx = 0, which is represented by the 
distribution 

d d d 



span 



dx ' ^ dz^ dy 

Lagrange-d'Alembert's principle gives the equations of motion 

X + yz — 
y = 
z — yx — 0. 

Discretize the system by defining the discrete Lagrangian Ld 

1 
2 



as 



Ldixo,yo,zo,xi,yi,zi) 



Xi - Xq 



yi - yo 



h 



zi - zo 



Then the discrete nonholonomic equations are 



X2 — 2X1 + 2^0 



Z2 - 2zi + Zq 

y2 - 2yi + yo 



Z2 - Zq 



X2 - Xo 



yi I ) = U (14a) 

(14b) 

Regarding as a Lie group under translations, the Euclidean metric is bi-invariant. 
Since L is induced by this metric and Ld is left-invariant, we have preservation of 
energy by Theorem [T] Figure [2] compares the energy behavior for our method 
against the DLA (discrete Lagrange-d'Alembert) algorithm in [TU] , 

In order to write the discrete nonholonomic momentum equation in Theorem |3] 
with respect to this group action, take two linearly independent sections of given 
by ^i{x, y, z) — (1,0, y) and ^2(2;, y, z) — (0, 1, 0). The equation for ^1 reads 



X2 - Xi 



2/2- 



Z 2 - Zi 



Xl - Xq 



yi- 



Z\ — ^0 
/l2 



(2/2 - 2/1) 



Z2 - Zi 



which turns out to be (14a I. Similarly, if we consider ^2 we reobtain (14b I. 

The DLA method proposed in [T^ also yields equations (14a I and (14b I, which 
is reasonable since both methods fulfill the discrete nonholonomic momentum equa- 
tion. However, the DLA method replaces ( 14c I by a discretization of the constraints 
that does not involve (xq, ?/0: ^0), such as 



Z2 - Zl 



2/2 + 2/1 \ X2- Xi 



= 0. 



Example 2. The snakeboard is a modified version of the traditional skateboard, 
where the rider uses his own momentum, coupled with the constraints, to move 
the system. The configuration manifold is Q = SE(2) x with coordinates 
{x,y,9,^p,(j}) as in figure [sj The center of the board, which is also the center 
of mass, is located at {x,y). We are considering here the case where the angles 
of the front and rear wheel axles are equal and opposite, as in [HJ However, 
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0.95 



Energy 



Our method 



DLA 
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50 60 



70 



90 100 



Figure 2. Energy behaviour for the nonholonomic particle using 
our method and the DLA method in pUj . 



we measure these angles with respect to the board instead of the x-axis. Figure |3] 
shows a configuration with all the angles positive. 




Figure 3. The snakeboard. The dashed line is aligned with the 
X-axis (not depicted). 

The continuous system is described by the Lagrangian 

L{q, q) = im(i2 + f) + i{J + 2J,)e^ + ^ Jo^^ ^ j^^2 

where m is the total mass of the system, J is the moment of inertia of the board 
about its center, Jo is the moment of inertia of the rotor mounted on the board 
and Ji is the moment of inertia of each wheel axle about its center. We assume the 
moments of inertia of the axles about the center of the board to be included in J. 
The distance between the center of the board and the wheels is denoted by r. 
The wheels are not allowed to slide sideways, so the constraints turn out to be 

X sm{9 + (f)) — y cos{9 + (fi) + rO cos{(j)) = 

isin(6' — 0) — ycos{9 — (/)) — r9 cos{(/)) — 0. 

If we define the functions a = rcos0cos0, b — rsm9cos4> and c — — sin0, then 
the constraint distribution is 

^ ( d d d ^ d d 

^ \ dip ^ d4)^ dx dy d9 

Endow Q with the Riemannian metric associated to the Lagrangian. This is repre- 
sented in coordinates by the diagonal matrix 

I = diag(rn, m, J', Jq, 2 Ji), 
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where J' ~ J + 2 Ji. The orthogonal complement to V is then 

^1 ( rf d d , d d 

V = span < J C7- ma — ,b— a- 

ox 00 



The projection Q : TQ 



dx dy _ 
T)^ is given in coordinates by the matrix 
J'c^ + mlP' —mab —J'ac 



-mab 
-mac 







J'c^ + ma? 
—mbc 





-J'bc 






which depends on (0, <j)), and its dual Q* is represented by the transpose. 

Consider the discretization of this system determined by the discrete Lagrangian 



LdiQk, Qk+l) 



1 
1 

2^ 



m{Axl 



Ayl) + l{J+2J,)A9l 



JiA( 



where qk = {xk,yk,Ok,i^kAk) (a column vector) and Azu = zu+i - Zk- 
The discrete nonholonomic equations Q can be written as 

D,Ld{qk,qk+l) + (Id - 2Q)l^{D2Ld{qk^i,qk)) = 0, 
SO in matricial form we get 

^ {lAqk + (Id - 2Ql){-lAqk-i)) = 0, (15) 

that is, 

qk+i = (Id - 2I-^Q^J)Aqk-i + qk- 

Regarding the configuration space SE(2) x as a Lie group, L is left-invariant. 
However, it cannot be right-invariant, because there are no bi-invariant metrics in 
SE(2). If one changes the group structure for the variables {x,y,6) from SE(2) to 
X S^, then both the continuous and discrete Lagrangians are bi-invariant. The 
numerical method itself does not depend on which symmetry group one takes, but 
considering this last group structure allows us to apply Theorem [T| to show that 
there is preservation of energy. 

On the other hand, we can still use the non-abelian group structure to write the 
discrete nonholonomic momentum equations, since only the left-invariance of Ld is 
required. Let us consider the action of the subgroup SE(2) on SE(2) x T^, and take 



62 



/o 
001 
Vo 



and 63 



/o -1 
100 
Vo 



Consider 



the typical basis of se(2): ei = || || 

the section |: Q ^ se(2) defined by |(x, y, 6, ip, 4>) = {a{9, (l))+c{9, (l))y)ei + {b{6, 0) — 
c{0, (j>)x)e2 + c{e, (/))e3, so we have |q ^ a{0 ,<!>)-§- + b{0, <!>)-§- + c{0, (p) ^ . Therefore, 



the discrete nonholonomic momentum equation in this case is 
(jf)|(9.,'Zfc+i)-(^f)|(9fc-i,9fc)- 

Xk+l - Xk 



m{a{9k+i,4'k+i) - a{0k,4>k)) 
+m{b{ek+i,(f>k+i) - b{0k, 0fe)) 



Vk+i - yk 



+ (J + 2Ji)(c(0fc+i,0fc+i) 



/l2 

c{0k,cj)k))- 



'fe+1 



/l2 



As an additional application, our method is ready to introduce controlled external 
forces. For instance we have added two controls: one applying equal but opposite 
torques on the wheel axles, and the other one on the rider. This was done by 
including appropriate terms on the right-hand side of equation (151. The figure 
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below shows a simulation where the snakeboard starts from rest and the controls 
are sinusoidal, with the same phase and frequency. This achieves the typical "snake- 
like" forward motion of the snakeboard, with increasing speed. 




Figure 4. The controlled snakeboard, moving left to right. 



Example 3. The Chaplygin sleigh consists in a rigid body that moves on a plane 
and is supported at three points. One of them is a knife edge and cannot slide 
sideways, and the other two can slide freely. Assume that the sleigh is symmetric, 
meaning that the center of mass is located on the line determined by the knife edge, 
at a distance a of the point of contact {x,y) (see figure [s]). 



Figure 5. The Chaplygin sleigh. 



The position of the sleigh is determined hy q ~ {x,y,9) E R'^ x S^, and the 
nonholonomic constraint is isin^ — ycosO = 0. If ni is the mass of the sleigh, / is 
its moment of inertia and {xc,yc) denotes the position of the center of mass, then 
the Lagrangian is 



L ^ (i^ + yc) + ^I9'^ = \ x- 



2a9x sin 



y 



The kinetic energy metric is represented by the matrix 



m 



m 




—am smu 
am cos 9 

—am sin am cos 9 I + ma^ 
so the constraint distribution and its orthogonal complement are 

d . „d d 



V — span < cos 



dx 



sm ( 



span 



sin 9 



dx 



dy' 89 

^Qd_ 

dy I + ma^ 89 



d 



The dual of the projector onto is then given by 



Q* = 



sm a 

— sin 9 cos ( 




' sm V cos ( 

cos^e 




/ + ma^ 
am cos 9 

I + ma^ 




2q2 



2a0y cos9 + a^9 



Discretize the Lagrangian by replacing x by (xi — xa)/h (analogously for y and 
, and 9 by (^o + Si)/2. We have applied the DLA algorithm, discrctizing the 
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constraints by {x2 — xi) sin((0i + 02)/2) — (y2 — J/i) cos((0i + 02)/2) = 0, and com- 
pared the results with the trajectory of the continuous system. This trajectory was 
obtained by applying standard numerical methods to the Lagrange-d'Alembert 
differential equations (see for example [2, p. 25]). Figure [6] shows the evolution of 
{{xk — + {Vk — VkY + (^fc — SkYY^"^ for both DLA and our method, where 
(Sfci yk^ Sk) are the values at t = hk of the trajectory of the continuous system. The 
results shown correspond to a particular trajectory with the initial points extracted 
from the continuous solution, but in general the errors are similar for the two meth- 
ods. We used m = J = 1, a = .2, 90 = (0,0,0) and qi = (-.2395, -.0070, .0589), 
which produces the heart-shaped loop typically described by the sleigh. 

It is worth mentioning that if we take a different discretization of the constraints 
for the DLA algorithm, such as {x2 — a;i)sin0i — (y2 — yi)cos(?i — 0, the error 
becomes larger by one to two orders of magnitude. Taking the right discretization 
is crucial in the DLA algorithm; in contrast, the accuracy of our method is close to 
that of DLA without the need of such a choice. 



7. Conclusions and future work 

In this paper, we propose a geometric integrator for nonholonomic mechanical 
systems for which the constraints are not required to be discretized. The integrator 
is different from the usual discrete analogue of the Lagrange-d'Alembert (DLA) 
principle which is presented in the works [TD] HI] • As initial conditions we propose 



points {qo,qi) satisfying (12 1 



Our method preserves in average the nonholonomic constraints, and the non- 
holonomic momentum map is also preserved. In addition, when the configuration 
space is a Lie group and some invariance conditions for the continuous and discrete 
Lagrangians are satisfied, we prove that the energy is preserved. In the particular 
case of a typical symmetric discretization of a mechanical Lagrangian we obtain 
a natural generalization of the well-known RATTLE method for holonomic con- 
straints. In addition, several interesting concrete examples illustrate these results. 

Of course, much work remains to be done to clarify the nature of discrete non- 
holonomic mechanics. A large part of this future work was stated in [52] and, 
in particular, we emphasize the following important topics: a complete backward 
error analysis and the construction of a discrete exact model for a continuous non- 
holonomic system; studying discrete nonholonomic systems that preserve a volume 
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form on the constraint surface, mimicking the continuous case; analyzing the dis- 
crete Hamiltonian framework; and the construction of integrators depending on 
different discretizations. 

For the case of reduced systems, it is possible to adapt the Lie-groupoid tech- 
niques introduced in the papers |16l 119] . considering now a fibred metric on the 
associated Lie algebroid and the induced orthogonal projectors. 

In future works, we will study these problems and, moreover, we will develop 
explicit constructions of higher order nonholonomic methods and applications to 
numerical methods for optimal control problems (of nonholonomic systems) . 
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